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Abstract. We consider the links between nonlinear dynamics and thermodynamics in the framework of a 
simple nonlinear model for DNA. Two analyses of the phase transition, either with the transfer integral 
approach or by considering the instability of a nonlinear particular solution, are discussed. Conversely, the 
computation of the largest Lyapunov exponent is obtained within a thermodynamic treatment. Differences 
with the Peyrard-Bishop model are also discussed. 
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1 Introduction 

The main goal of the paper is to study the links between 
the microscopic dynamics and macroscopic thermodynam- 
ical properties in a very simplified model for DNA. Both 
aspects are usually not studied simultaneously; in the lit- 
erature, the main goal is often, either to consider dynami- 
cal aspects of coherent structures (solitons for example) in 
a system at zero temperature, or to derive thermodynami- 
cal properties without really considering the consequences 
of the existence of these coherent structures. Here, on the 
contrary, we will put the emphasis on the link and show 
that both approaches give important insights to the de- 
scription of the physical properties. 

Two aspects will be of particular importance and we 
would like to shed light on them already in the introduc- 
tion. A recent method [Hj, which showed how the stabil- 
ity of a nonlinear solution of the dynamical equations ex- 
hibits an interesting approach to describe a phase tran- 
sition, will be applied in this new model: it will reveal 
how the small amplitude dynamics needs a careful treat- 
ment to adequately describe the thermodynamics. Second, 
we will describe the relationships between the maximum 
Lyapunov exponent which characterizes the dynamics and 
the thermodynamics. Using a geometric method [21], we 
will explicitly compute the evolution of this dynamical 
quantity as a function of the temperature: consequences of 
the phase transition will therefore be explicit. From these 
calculations, we predict features absent in the Peyrard- 
Bishop model, although both models are very similar: a 
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non monotonic behavior of the maximum Lyapunov ex- 
ponent in the low temperature phase, and a jump at the 
critical temperature. 

In Sec. |3 we present briefly the model. In the following 
Sec. we show how one can derive its thermodynamical 
properties. Then, the emphasis is put in Sec.^to a special 
particular solution and on its use to explain the thermody- 
namical properties. Finally, using the powerful geometric 
method introduced recently to compute the largest Lya- 
punov exponent, we discuss in Sec. El the link between this 
dynamical quantity and the phase transition, a thermody- 
namic concept. 



2 The Model 

A very simplified model has been proposed in 1989 by M. 
Peyrard and A. R. Bishop PP to describe DNA denatura- 
tion. This biological phenomenon leads to the breaking of 
the H-bonds linking both strands of DNA. This process, 
which appears when either the temperature increases or 
when the pH of the surrounding solvent is modified, was 
previously successfully described by Ising models • How- 
ever, the dynamical properties of the phenomenon were 
not captured by these static descriptions. This is why, 
keeping the principle of simplicity, Peyrard and Bishop de- 
scribed this highly complicated biomolecule by two chains 
of particles coupled by nonlinear springs. This step to- 
ward a more complex system, since including the minimal 
dynamics, was unexpectedly more successful than first un- 
derstood and have lead to many studies (See Ref. [3] for a 
review). A large part of the results were interesting from 
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the physical point of view: one may list in particular stud- 
ies of discrete breathers modes and energy localization 
in systems involving nonlinear and discreteness effects fjl 
E], existence of phase transition in one-dimensional sys- 
tem . . Furthermore, several recent results have empha- 
sized that this model could be successfully used to locate 
promotor regions jZj for real DNA sequences. This un- 
expected report and similar ones explain why so simple 
models are still nowadays thought to be possible powerful 
tools to describe real DNA dynamics. 

Two linear chains of particles describing phenomeno- 
logically the nucleotides describe the two different strands 
of DNA as schematically represented in Fig. ^ In the sim- 
plest description, all particles of the same chain are har- 
monically coupled whereas the interstrand interactions are 
restricted to facing nucleotides; long-range interactions are 
neglected at this level of description. It is important to un- 
derstand that because the goal is to describe denaturation 
dynamics, only transversal degrees of freedom are taken 
into account. Defined with respect to their equilibrium po- 
sition, the displacement of the center of mass of the n th 
nucleotide are called u n (resp. v n ) for the top (resp. bot- 
tom) chain. 
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Fig. 1. Schematic presentation of the DNA model. 

Denoting p u the conjugated momentum to the spatial 
position u and the number of nucleotides being TV, the 
Hamiltonian model can be written as 
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The dynamics and the thermodynamics of the first 
part, H x , which corresponds to a linear chain of harmonic 
oscillators, can be easily computed. We will omit this part 
in the remaining of the paper without loss of generality. 
The second part, H yi on the contrary needs further de- 
velopments. For the sake of simplicity, we will omit the y 
index. 

The system defined by this hamiltonian H exhibits a 
second order phase transition as shown in the next sec- 
tion. The low temperature phase corresponds to states 
where the particles are located close to their equilibrium 
position, the associated DNA being in the native state: it 
must correspond to the bottom well of the potential V. On 
the contrary, in the high temperature phase, DNA being 
denaturated, the link between facing nucleotides of both 
strands are broken: consequently, associated particles of 
the model must be located in a plateau of the potential, 
far from their equilibrium positions since the force is van- 
ishing. The position in this plateau will be thermodynam- 
ically chosen because of the entropy contribution to the 
free energy, important only at high temperature. 

This simple physical description guides therefore the 
appropriate choice for the potential V. Nevertheless, the 
analytical calculations that we will present now are pos- 
sible for only a few cases. The Morse potential V m (y) — 
£>(exp(— a m y) — l) 2 was the first choice 0E1- Here, we will 
present another possible example [B]: 



2 — if V > 

V(y) = i cosh ay 

+oo if y < 0. 



(4) 



As we will see, the qualitative shapes are really close but 
the differences of curvatures lead to several consequences. 
In addition, the impossibility to reach negative values for 
the variable y gives interesting properties. One important 
question concerns the influence of the details of the poten- 
tial V. We will discuss in particular the following points: 

i) values of the critical temperature of the phase transi- 
tion for both potentials, 

ii) consequences to the related characteristic lengths: or- 
der parameter and correlation length, 

iii) largest Lyapunov exponent as a function of tempera- 
ture. 



The interstrand potential V describes the effective interac- 
tions, i.e. in particular hydrogen bonds between base pairs 
but also the repulsion between phosphates. The canonical 
transformation x n = (u n + v n ) / s/2 and y n = (u n ~ v n )/^/2 
decouples simply both degrees of freedom since H can be 
rewritten as H = H x 



3 Thermodynamics of the model 

3.1 The canonical partition function 

As it is well-known nowadays the statistical me- 

chanics of such a one-dimensional short-range Hamilto- 
nian can be exactly derived with the transfer integral 
method (See appendix 0). In the framework of the con- 
tinuum approximation, the solution relies on solving the 
following Schrodingcr equation 
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Fig. 2. Comparison between potential Q (represented by 
the solid line) with appropriate characteristics (a, j/o) and the 
Morse potential V/D (dotted line). 



if the lattice spacing between sites in the x-direction is set 
to one. 

Defining the quantity 




where 



VKD 



(6) 



(7) 



this equation has N„ = E(jq + 1/2) localized states, with 
E(.) denoting the integer part. One notes that T c corre- 
sponds to the disappearance of the last discrete state, and 
will be called the critical temperature. 

The (N v + 1) localized states tpk can be expressed jS] 
in terms of hypcrgcometric functions as 



~>Pk{y) 



cosh(a?/) b2fc + 1 
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where A/fc is the normalization factor of the wave function, 
b n = 2r\ — n and finally 



Sk = —. 



2(3 2 K 



(2?7 — 2fc — 1) 



(9) 



the associated eigenvalues. 

The ground state which is particularly useful in the 
remaining of the paper can be simplified as 



■00 (y) = 4 /2a 
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(10) 



by introducing the Euler function 
r(x + y) 



B{x,y) 



r(x)r( y ) 



dt(l-t) *~H V 



(11) 



3.2 Discussion of the choice of the parameters set 

One of the goal of this work is to make a detailed compar- 
ison between the Morse potential and potential Q. The 
parameter set is obtained either by considering the effec- 
tive physical interactions, or by choosing the values to get 
the same melting temperature. Assuming that the depth 
of the potential is known, both possibilities will define a 
relation between the width of the potentials: for the 
Morse one, and a -1 for potential Q. 

The choice a m « 1.472 a corresponds to the best fit. 
On the contrary, as the transition temperature for the 
Morse potential flTTTll9] is T™ = T C V8, the second choice 
a m = ay/8 would lead to the same critical temperature. 
Consequently, when the parameters are fitted, the critical 
temperature differs by a factor two! 

This unexpected disagreement reveals a hidden differ- 
ence between both potentials: the underlying reason is 
the important contribution of negative positions y, in the 
Morse case. They cannot be neglected even if the fast ex- 
ponential increase was thought to play the role of the im- 
possibility of interpenetration. 

The appropriate solution is to introduce an additional 
parameter. The inverse width a being given by the critical 
temperature, one could adapt the minimum of the poten- 
tial to minimize the differences with the Morse potential. 
A fit restricted to positive j/-values over both variables 
(a m , yo) leads to an excellent agreement. This case is pre- 
sented in Fig. by the solid line. 

All results below correspond to the following set of pa- 
rameters to = 300 u. a., D = 0.00094 cV, K = 1.9 eV.A~ 2 
and a = 4.5 A . The value of T c will be different from 
the Morse case, but both potentials will be very similar 
as emphasized by Fig. [5] This choice will allow a precise 
study of the negative y-values region and of the impor- 
tance of potential curvature for the Lyapunov exponent 
discussed in section [5] 



3.3 Characteristic lengths 

As usual, the thermodynamic properties of this system can 
be characterized by an order parameter, and its fluctua- 
tions, in the vicinity of the critical temperature T c . Here, 
the appropriate choice is the quantity 



dx x | iI)q(x)Y 



(12) 



which diverges for T = T c . This clarifies the name critical 
temperature which separates, the phase with a finite order 
parameter (native state) from the phase with infinite order 
parameter, representing the denaturated state. 
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Fig. 3. Order parameter al (solid line and triangles) and cor- 
relation length a£± (dotted line and diamonds) as a function 
of the reduced temperature t = 1 — T/T c . Lines corresponds to 
the transfer integral result within the continuum approxima- 
tion (Eqs. (11 21 and dJ) whereas symbols correspond to the 
exact numerical solution of the operator 1571 . Presented in the 
inset, the exact numerical data plotted with logarithmic scales 
emphasize the critical behavior close to the critical tempera- 
ture. The dashed (resp. dash-dotted) line corresponds to the 
asymptotic expression (1131 (resp. (1151 ) . 



The associated critical exponent can even be deter- 
mined by using the asymptotic expression J Q dx x/cosh" x 

I /a 2 . Using expression (|10J) and introducing the usual 
reduced temperature t = 1 — T/T c , one finally gets 



3 1 1 
16a|t| 



(13) 



In agreement with critical phenomena theory we obtain 
therefore the usual critical exponent (3 — — 1, for this sec- 
ond order phase transition. Fig. |21 presents the evolution 
of £ as a function of the temperature. Formula 112f) . rep- 
resented with a solid line, agrees very well with the exact 
result obtained with the discrete transfer operator (|57jl . 
In the inset, the logarithmic plot emphasizes the critical 
behavior and confirms the scaling exponent f3 = — 1. 

Interestingly it is also possible to determine the fluc- 
tuations of the order parameter which can be computed 
from the following expression 



(14) 



As above, the asymptotic expression can be easily derived 



by taking into account that J °° da; x 2 /cosh" : 



2/a 3 



One thus obtains in the vicinity of the critical temperature 
that 

+ t c 3^3 1 1 



6 



(15) 



16 a \t\ 

The critical exponent is therefore v± = 1. Fig. [21 shows 
that the order parameter and its fluctuations are of the 



same order in this temperature regime. The inset confirms 
again the critical behavior and the scaling exponent v = 1. 

It is important to stress that all above analytical re- 
sults are confirmed by the numerical but exact solution 
of the transfer integral operator (|57|) . without relying on 
the continuum approximation: see symbols in Fig [3] As it 
is nowadays known [EK], the critical exponents are valid 
even in the region where discreteness effects can not be 
neglected. To be more specific, let us introduce the pa- 
rameter 

(16) 

measuring the onsite potential with respect to the elastic 
coupling. For small value of this parameter, the contin- 
uum approximation is valid and results (|13|) and (|15(1 con- 
firmed. For the set of parameters chosen in this work, R is 
equal to 10 -2 , and one get results in values slightly differ- 
ent from Eqs. I|13(l and (|15f) for T c , £, and However, as 
shown by Fig. [31 the differences are hardly distinguishable 
and moreover the critical exponents are fully confirmed. 



4 Domain Wall 

In the framework of this model, it is also interesting to 
discuss a new method jSJEl i alternative to the usual ther- 
modynamic one proposed in previous section, to detect 
the phase transition from dynamical considerations. This 
will illustrate again the importance of the wall located at 
y = 0. 

Equations of motion corresponding to Hamiltonian Q 

are 

dV 

my n = K(y n+1 + y n _ x - 2y n ) - - — (17) 



dy r , 



or, in the continuum limit, 



.. _ K d 2 y 1 dV 
m dx 2 m dy 



(18) 



The uniform profile at the minimum of V{y) is a static 
solution of the infinite chain with free ends. Profiles veri- 
fying d 2 y/dx 2 = are approximate solutions only on the 
plateau of the potential, since dV/dy is close to zero for 
large y. 

There exists also an exact, unbounded, domain-wall 
like solution 



- ^-Argshe ±z = — In 



Vdw( x ) 



3 ±z 



= ±2z 



(19) 



where z = \2R(x — xq) and xq is an arbitrary constant. 
Solution (|19J) represents a configuration which links the 
stable minimum to a particular member of the mctastablc 
configurations, with a slope ^J2D/K. One can easily checks 
that this corresponds to equal contributions to the elas- 
tic and the on-site potential energy densities (D per site). 
Consequently, the energy of the solution contains a term 
which is proportional to the number of sites to the right 
of xq and if lattice sites are numbered from to N, one 
has 



E 



DW 



(N - x )2D + O(N ) 



(20) 
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At zero temperature the profile 1)19(1 is consequently not 
stable, and the wall spontaneously move to the right, "zip- 
ping" back the unbound portion of the double chain. This 
instability changes however under the influence of temper- 
ature. 

At non zero temperatures, let us consider small devi- 
ations with respect to l|19|) . i.e. 



y(x,t) = y^ w (x - x Q ) + ajfj( x ~ xo)e 



(21) 



where \ctj\ -C a 1 . The linearized eigenfunctions fj satisfy 
the Schrodinger-like equation 



dz 2 



1 -2e 



2.: 



(1 + e 2 *) 



fj 



,2 J J 



2KR 



(22) 



Eq. 1)22(1 has no bound states [S]. There are however scat- 
tering states: acoustic phonons oscillating on the flat por- 
tion of the potential with frequencies 



IK . 

(1 — cos q) 

m 



(23) 



are some of them. 

In the bottom of the well of the potential, let us first 
forget the wall. Consequently scattering states would be 
optical phonons with frequencies 



J opt 



2D a 2 IK. 1 

i (1 - cosq). 



(24) 



At finite temperatures, the domain wall would there- 
fore be accompanied by a phonon cloud contributing to 
the free energy as 



F ph = k B Tx 



(25) 



where we omit terms independent of xq. Introducing the 
dispersion relations, we can evaluate |15j the integral in ((25(1 
using /J'dxliifl - cos a; + R] = 2?rln[(v r R+ VR + 2)/2]. 
We obtain thus the total free energy (DW plus phonon 
cloud) 



F = k R T\n 



VR+VRT2 



V2 



2D]x a + const. (26) 



This result describes in very simple terms why and when 
the phase transition occurs. At temperatures lower than 



T c = 



2D 



ks In 



^R/2+^l + R/2 



(27) 



the prefactor of xq in Eq. ((26(1 is negative, describing the 
DW's natural tendency towards high positive values of Xq: 
it "zips" the system back to the bound configuration. Con- 
versely, at temperatures higher than T c , thermal stability 
is achieved and the DW "opens up" . 



It should be noted that the value of T c predicted by 
the above DW argument coincides exactly with the result 
obtained for the Morse potential [HE]. The limiting be- 
havior in the continuum approximatio n, i.e. in the limit 
R < 1, leads for example to T c = 2\f2KD/(ak B ). This 
result differs nevertheless from Eq. Q. 

Expression ((27(1 is consequently not valid for the sech 2 - 
potential J3J of interest. The underlying reason is the im- 
possibility of usual phonons to take place in the bottom 
of the well. The harmonic approximation of the problem 
leads indeed to a nonlinear problem because of the non- 
linear condition y > introduced by the wall. Optical 
phonons with frequencies 1(24(1 are therefore totally modi- 
fied by the wall and cannot be computed. A way to take 
into account this nonlinear condition would consist in se- 
lecting only the even modes of optical phonons l(24|) . How- 
ever, it turns out to be unexact and unsufficient . This il- 
lustrates again that the presence of the wall strongly mod- 
ifies the dynamics of the system, which leads to important 
thermodynamic consequences. 



5 Geometrical method to derive the largest 
Lyapunov exponent 

In the theory of dynamical systems, the concept of Lya- 
punov exponent has also attracted a lot of attention |16l 
117) because it defines unambiguously a sufficient condi- 
tion for chaotic instability. Unfortunately, except for very 
few systems, it is already an extremely difficult task to 
derive analytically the expression of the largest one, Ai, 
as a function of the energy density. As some promising re- 
sults have been recently obtained to describe some prop- 
erties of high-dimensional dynamical systems 18.19,20. 
I21II22II23) , by combining tools developed in the framework 
of dynamical systems with concepts and methods of equi- 
librium statistical mechanics, the idea that both concepts 
could be related was proposed (24) . 



5.1 Riemannian geometry approach 

The main idea is that the chaotic hypothesis is at the ori- 
gin of the validity of equilibrium statistical physics, and 
this fact should be traced somehow in the dynamics and 
therefore in the largest Lyapunov exponent. The method 
is based on a reformulation of Hamiltonian dynamics in 
the language of Riemannian geometry (24) : the trajectories 
are seen as geodesies of a suitable Riemannian manifold. 
The chaotic properties of the dynamics are then directly 
related to the curvature of the manifold and its fluctu- 
ations. Indeed, negative curvatures tend to separate ini- 
tially close geodesies, and thus imply a positive Lyapunov 
exponent; nevertheless chaos may also be induced by pos- 
itive curvatures, provided they are fluctuating, through a 
parametric-like instability. To approximate the curvature 
felt along a geodesic, the method uses a Gaussian statisti- 
cal process. The mean value of this process is given by kq 
and its variance by a K , where kq and cr K are the statistical 
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average of the curvature and its fluctuations, which can be 
computed by standard methods of statistical mechanics. 

Finally one ends up with the following expression of 
the largest Lyapunov exponent [21] 



A - 1 ( A 4K ° 



where 



A = 




(28) 



(29) 



In this definition, r, the relevant time scale associated 
to the stochastic process, is function of the two follow- 
ing timescales: t\ ~ 7r + a K is the time needed to 
cover the distance between two successive conjugate points 
along the geodesies, whereas r 2 ~ t/koI<J k is related to 
the local curvature fluctuations. The general rough phys- 
ical estimate r ~ (1/ti + l^)" 1 completes finally the 
analytical estimate of Ai. We continue now by calculating 
the mean value of the curvature and its fluctuations as a 
function of the energy density. 



5.2 Average curvature 

The curvature of the Ricmannian manifold is directly given 
by the Laplacian of the potential. One needs therefore to 
compute the microcanonical average of the quantity 



AV = 2KN + 2a 2 Dj2d(yk) 



(30) 



and the mathematical formula I(a+l) = I(a)2a/(2a + 1). 
canonical averages can be simplified. Introducing the pa- 
rameter R defined in Eq. (|16(l . one gets 



K 



N 



2K [1 + 4R 



(277 - 1X277 - 3) 
(47?+1)(4t7 + 3) 



(35) 



Above expression gives in particular the following limiting 
behaviors 

2Da 2 (36) 



lim k = 2K 



and 



lim «n = 2K , 



(37) 



which coincide with asymptotic results for the Morse po- 
tential [H| . 

It is however important to notice that expression (|35J) 
suggests that the mean curvature is not always positive, 
its sign being tuned by the value of R. As negative curva- 
tures enhance dynamical instability, this result may have 
strong consequences on the largest Lyapunov exponent. 
A careful study jT2] , shows that for values of R larger 
than R c = (31 + 3\/105)/8, expression lT3*5)l could be neg- 
ative in a given interval of temperatures. This result has 
to be criticized since expression <|35|) has been derived in 
the continuum approximation, and one expects important 
discreteness effects for R values as large as R c . Numeri- 
cal, but exact, resolution of the transfer integral operator 
shows that the curvature is actually positive for all R. 
Fig. pj] shows the curvature as a function of temperature 
for the parameter set chosen in this work. 



5.3 Fluctuations of the curvature 



and its corresponding fluctuations, where 

q i 

g(y) 



cosh 4 ay cosh 2 ay 



We finally obtain 

{AV)p = N (2K + 2a 2 D{g{y k )) ll ) 



(31) 



(32) 



As, in the thermodynamic limit, ensemble equivalence 
ensures that averages are equal in all statistical ensem- 
bles, we will compute them in the canonical one since the 
transfer operator method has been shown to be a power- 
ful tool to compute thermodynamic functions, especially 
if the continuum approximation is valid. 

Calculation of Eq. 132f) relies on the computation of 
(g{ay)) ca _ n , i.e. of terms such as (1/cosh a ay) can . Using 
the transfer operator method in the continuum framework, 
we immediately find 



1 



cosh 2 " ay 



dy X \Mv)\ 2 - (33) 

o cosh ay 



With expression ^Bj 
1(a) 



o cosh a x 



da; „ /„ 1 



(34) 



Contrary to the statistical averages such as the curva- 
ture Ko, fluctuations are ensembles dependent. This is why 
one computes first fluctuations in the canonical ensemble, 
before using the Lcbowitz-Percus-Verlet formula J3| to 
get the microcanonical fluctuations. 

In the canonical ensemble, using Eq. (|61|l . fluctuations 
of expression (|30() are given by 



(S 2 AV) 1 



Aa 4 D 2 N N 



Y,^(y l )g(y ] ))~(g(y l ))(g(y 3 ))} (38) 



^(g(yN)g(yN-i)) - N(g(yY 



(39) 



EE' 

I q=0 



-0l(e g -e o ) 



= e -0l(e q -e o ) 
I q=l 



N^oS— 1 1 - e -0(e q -ea) 
9=1 



dyg(yW q {y)ipo(y) 
-N(g(y) 2 ) 
dyg{yW q {y)i>o{y) 

dy g(y)i>a(y)^o(y) 



(40) 

I 

'(41) 



(42) 



Above expression can be used to compute the canonical 
fluctuations. 
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Fig. 4. Average curvature kq and canonical fluctuations a K 
versus temperature. The triangles (resp. diamonds) represent 
the exact numerical result for the curvature (resp. canoni- 
cal fluctuations), obtained with the transfer integral formal- 
ism. The solid line corresponds to analytical formula il'-iol 
whereas the dashed line correspond to the microcanonical fluc- 
tuations ijl^. 



The microcanonical fluctuations will finally be recov- 
ered by using the Lebovitz-Percus-Verlet formula For 
any quantity C with fluctuations (S 2 C), both fluctuations 
are related through the formula 



<5 2 C*) M = (<5 2 C) can + 



d(U) a 
8(3 



1 f d(C) 



03 



(43) 



where (U) ca n is the averaged energy of the system. 

The canonical partition function Z being the product 
of a kinetic part Zt and a configurational one Z c , the 
averaged energy is given by 



<£/>c 



din Z T dlnZ c 



dp 



0ft 



(44) 



The first contribution is as usual N/ {2(3) whereas the last 
one, can be simplified in the continuum approximation by 
using the transfer integral method (See Appendix lAl and 
in particular Eq. I|fi0[l). Denoting £q the ground state of 
the associated Schrodingcr equation, one obtains In Z c = 
-N/3e with 



1 ,„ @K a 2 2 



The averaged energy is therefore 



(C/)ca„ = N 



1 



d(Pe ) 
00 



(45) 



(46) 



The final analytical expression for the microcanonical fluc- 
tuations is not simple. Introducing quantities (3 C = I/QcbTc) 
and S = T c /T, one gets 



(<5 2 «o) M - (S 2 k ) c 



l8432a 4 D 2 K(3 c 



= -S 5 (-165 2 + 7y/l + 8S 2 



/(l + 8<5 2 ) z/[2 + Vl + 85 2 
/(16K(3 c S 3 y/l + 85 2 - 36a 2 5 2 
+40a 2 S 2 Vl + 8<5 2 + 2K(3 C 6 
3a 2 



Vl + 85 2 



5a 2 \/l + 8(5 2 ).(47) 



Above expressions (|42f) and l|47|) can be combined to com- 
pute the microcanonical fluctuations. This is what has 
been performed to plot the fluctuations of the curvature 
in Fig. 2| 

Close to the critical temperature T c , the continuum 
part of the transfer operator spectrum should be taken 
into account and an explicit analytical calculation is pos- 
sible in principle but particularly tedious. On the contrary, 
one can simplify above expression in the low temperature 
regime as shown in next section. 



5.4 Low temperature estimate 

In the low temperature regime by replacing the prefactor 
(1 — exp[— /3(s q — £o)]) _1 by 1 in Eq. l|4"2|l. one finally gets 



(S 2 AV) can ~ Wa 4 D 2 J2 (^Mi>o) (M9\%) 

9=1 



4Na 4 D 2 



(4>o\9 2 \1>o) - Wo\g\il>o) 



Combining result 



(V>o|3 2 |^o) = 16- 



277(277 — 1)(4?7 2 - 677 + H) 



(477 +l)(477 + 3)(4?7 + 5)(477 + 7) 
with formula l|35|l . one ends up with 

(277 - l)(256/7 3 - 184?7 + 105) 



48 



(4?7 + l) 2 (4r; + 3) 2 (4?7 + 5)(4r7 + 7) ' 



(48) 
.(49) 

(50) 
(51) 



As the microcanonical correction l|47|) can be neglected 
in this region, the microcanonical fluctuations of the cur- 
vature are finally given in the low temperature regime by 



(S 2 AV) C 



N 



AD 2 a 4 



48(27? - 1)(256t7 3 - 184t7 + 105) 
(4t; + 1) 2 (477 + 3) 2 (4?7 + 5) (4?/ + 7) 



(52) 
(53) 



This expression can be used in the low temperature re- 
gion to get a simpler expression for the largest Lyapunov 
exponent. One thus obtain that it increases quadratically 
with the temperature. 
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5.5 Largest Lyapunov exponent 



6 Conclusion 



We can now estimate the largest Lyapunov exponent Ai 
for this high-dimensional system of N coupled particles in 
the external sech 2 -potential J3J. 

As kq is positive, the instability of trajectories is due to 
fluctuations of curvatures as reported by Pettini, Casetti 
and Cohen [21] . Analytical expressions (|35|l for kq and l|53|l 
for a K were the only missing points in the estimate 1|28[) for 
A\. One can in addition notes that in the low temperature 
region, where <j K <C k as attested by Fig. lH one gets the 
asymptotic expression 



Ai 



K - 



(54) 



We have presented a new qualitative model for DNA de- 
naturation directly inspired by previous works [T|H l|f(T| . Its 
complete statistical mechanics was derived, as well as all 
features related to the second order phase transition: not 
only the critical temperature, but also the critical expo- 
nents related to the order parameter and the transversal 
correlation length. 

We have in particular emphasized the important role 
of the negative y- values for the Morse potential which were 
believed to be unimportant fTI] ■ If critical exponents are 
of course not affected, the critical temperature is strongly 
dependent on it. Furthermore, using a geometric approach 
to estimate the largest Lyapunov exponent, we have com- 
puted its evolution as a function of the temperature. The 
results are unexpectedly qualitatively different from those 
obtained with the Morse potential. 



6x10 




0.4 0.6 
T/T c 



Fig. 5. Largest Lyapunov exponent Ai versus T/T c . The tri- 
angles correspond to the transfer integral result, whereas the 
dotted line correspond to the analytical expression 1541 . 
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A Transfer integral method for the canonical 
partition function 

With periodic boundary conditions, yo — yjv, the configura- 
tional partition function of Hamiltonian 10 can be written as 



Z c = 



by introducing the symmetric function 

1 



f{.Vn,yu-i) 



D 
~2 



<%o - Vn) 



1 



(55) 



cosh ay n cosh ay n -\ 



K , ,2 
+ y(?M ~ Vn-l) ■ 



Defining the transfer operator T as 

T[0](v)= [ dxcf>(x)e- 0fM , 



(56) 



(57) 



Fig. (JSJ) presents the temperature evolution of the largest 
Lyapunov exponents. Two features should be emphasized 
in comparison to what has been previously reported for 
the Morse potential f^l] . One can notice a local maximum 
and a local minimum of the Lyapunov exponent in the low 
temperature region. More importantly, one has to realize 
that for temperatures larger than the critical one T c , par- 
ticles are on the plateau of the potential; the chain being 
equivalent to a linear chain, the largest Lyapunov expo- 
nent Ai if of course zero. The Lyapunov exponent should 
thus present a jump close to the critical temperature. It is 
important to check these two strong predictions by consid- 
ering careful microcanonical numerical simulations with 
large systems. This would provide a precise test of the 
geometrical method to calculate Lyapunov exponents. 



its eigenvalues Ek and normalized eigenvectors ipk are given by 
T[ipk] = exp(-/3e fc ) ip k . 

Introducing the orthonormalization condition 



% - yo) =^2^Uyo)ipk(v) 

k 

in Eq. I|55|l . one gets |^] 



-N fit 



(58) 



(59) 



If the lowest eigenvalue Eo is discrete and located in the gap 
below the continuum, one can simplify above expression in the 
limit N — > oo, so that 



Z c ~ e 



-N0e o 



(60) 
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A similar method to compute the canonical average of any 
function h(y) leads to the following result 



{h{VN)c 



dyh(y)\My)\ 2 - (61) 



Above results are valid without approximations. However, 
applying the continuum approximation, it is possible to go one 
step further since there is a mapping |1 1191 between the transfer 
integral operator and the following Schrodinger equation 



l d 2 v> 



D 



2f3 2 K dy 2 cosh 2 ay 



tp = Skip- 



(62) 



As the spectrum of a quantum particle in the potential Q is 
known, one can derive the analytical expression of 161H within 
this approximation. 
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